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Abstract 


In this paper the numerical robustness of four generally applicable, recursive, least- 
squares- estimation schemes is analyzed by means of a theoretical round-off propagation 
study. This study highlights a number of practical, interesting insights of widely used 
recursive least-squares schemes. These insights have been confirmed in an experimental 
verification study as well. 
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1 Introduction 

Recursive least squares (RLS) schemes are used in a broad class of practical applications, 
such as in the self-tuning regulator schemes developed by Astrom et al. (1977). The 
simplicity of this class of schemes has been an important motive for their use. Some of the 
other reasons are 1) the versatility of implementations, allowing to exploit, special (shift) 
structures in the regression model, first demonstrated by Levinson ( 1947), 2) the well-known 
statistical properties of the RLS estimates, and 3) robust ways to cope with time variations 
in the regression parameters, see T. Hagglund (1983) or R. Kulhavy (1985). 

However, the numerical robustness of a number of RLS implementations is still not 
well understood. This can be concluded from the appearance of simulation studies reporting 
the divergence phenomenon in the implementation referred to in this paper as the conven- 
tional RLS implementation. Here the divergence phenomenon is a generic term used to 
indicate a whole class of problems where the quantities updated in the actual computer 
implementation lose their theoretical properties. A well-known example is the loss of sym- 
metry of the parameter error covariance matrix. 

Until now no explanation has been given for this phenomenon, only a number of 
“cures” have been proposed. These ernes have two disadvantages: 1) they might unnec- 
essarily complicate the numerical implementation and 2) they might give “satisfactory” 
results in simulation, but actually fail in real-time operation. 

In numerical analysis, an error analysis is performed to understand the numerical 
robustness of algorithmic implementations. For RLS, such an analysis has been performed 
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in S. Ljung and L. Ljung (1985). Although, this analysis reveals some fundamental un- 
derstanding of the robustness of RLS schemes, it does not answer the following important 
practical questions: 

1. When does the loss of symmetry occurs in the conventional RLS scheme? 

2. How does the loss of symmetry interferes with the loss of positive definiteness of the 
parameter error covariance matrix? 

3. Does the conventional RLS scheme result in numerically less-accurate results com- 
pared to square root type of RLS? 

4. When using a square root implementation do we choose a covariance- type or an 
information-type RLS? 

To answer these questions, a new and detailed error analysis is performed in this 
paper. This analysis considers the following four generally applicable implementations of 
RLS: 

1. Two implementations of the conventional RLS (indicated respectively as the CLSl 
and CLS2). 

2. The square-root-covariance RLS implementation (SC’LS). 

3. The square-root-information RLS implementation (SILS). 

These implementations are precisely defined in section 2. Section 3 then presents the results 
of the round-off error propagation study for these four implementations. The new insights 
from this error analysis are evaluated experimentally in section 4. Finally, section 5 presents 
some concluding remarks. 
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2 Recursive Least Squares Estimation 


2.1 The Linear Least Squares Problems 

Let the linear scalar regression model be denoted as: 

y k = + ek ( 1 ) 

where the regressor vector ^ € R n , e* is a zero mean discrete white noise sequence with 
variance <t\ 2 and (.)' denotes the transpose. When the observations of (yki'pk) have been 
obtained for k = 1 , • • • , JV (with N > n), the least squares (LS) estimate of d is defined as: 

N 

$N = arg ro*n]TA iV ~‘(y, - dV») 2 ( 2 ) 

1=1 

Here {A jV_I } is a sequence of weighting coefficients that result if we discount old data by 
so-called ‘‘exponential forgetting (A < 1).” For this case, simple calculations show that the 
LS estimate (equaion 2) is given by: 

#N = Rn 1 /n ( 3 ) 

R N = YlWi X ‘ W " ( 4 ) 

1=1 

s 

= (5) 

i=l 

For real-time operation, it is of primary importance to compute the estimates recur- 
sively. Therefore, the set of equations (3-5) is rearranged by simple manipulations to the 
following recursive form: 

$Jlr = 1 + R k l 'rk(yk - l^.-1'rfc) (6) 

R k = \Rk-\ + •Pk'rh (~) 

The set of equations (6-7) is updated from A- > 0 on. This can be done when for k = 0 the 
initial conditions «?_i and /?_i are specified. 

"In practice a is generally not known a priori, but is estimated during the operation of the RLS. However^ 
what actual value for <Ji is used in the RLS will turn out to be not important in the analysis of the numerical 
robustness 




2.2 The Conventional RLS Method (CLS) 


The recursive relationships (equations 6-7) are not well suited for computation, since the 
nxn information matrix R k has to be inverted each time step. Therefore, it is more natural 
to introduce 

p k = r; 1 


and by applying the matrix inversion lemma. Ljung and Soderstrom ( 1983), (equations 6-7) 


become 

dfc = + Kk{yk - 

Kk = 

\<r\ + ip^Pk-nPk 

P _ 1 ( P _ Pk-lfkfkPlr^l \ 

fc A \ A <rl + APu-i*J 


( 8 ) 

(9) 

( 10 ) 


Remark 1: The linear regression model (equation 1) can be written in the state space 
form, 

(ID 

yh = + e* with E (t*) = <rl (12) 

Together with the initial conditions (ii_i,P_i), this model is a special case of the general 
state space model used in the Kalman filter design (as described in Verhaegen and Van 
Dooren, 1986): 


fjk+i = A k *k + B k w k with £(*<>) = - E(x 0 x' 0 ) = P. x and E(w k w' k ) = Q k (13) 

y k = C k x k + v k with E(v k v k ) = R k (14) 

for the following particular system matrices: 


A k = I . B k — Q k — 0 , C k - <r' k and R k - <r k 


(15) 


For these system matrices (equations 8-10) represent the so-called measurement update of 
the conventional Kalman filter as defined in Verhaegen and Van Dooren (1986). Therefore, 
equations (8-10) are conformally indicated in this paper as the conventional RLS (CLS). □ 

A conunonly used implementation of the CLS is represented in Table 1 as the CLSl. 
Here we observe that for the calculation of the matrix K k the symmetry of P k - i is exploited. 
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The only reason for that is a slight reduction in computational complexity. Furthermore, the 
performed operation to calculate Pj‘_ l is inherently symmetric. Therefore, this implemen- 
tation “seemingly” forces the parameter error covariance matrix Pk to remain symmetric. 
However, the theoretical error analysis will prove that exactly the contrary is the case. 
Namely, that when implementing equations (8-10) exactly as they are, as done in the CLS‘2 
implementation of Table 1, the sensitivity to the loss of symmetry of Pk vanishes. 


quantity 

mathematical expression 

CLSl 

CLS2 

fiP 

SkPk - 1 

'f'k * p k- 1 


r i 

X&1 + <fi'i t Pk-i i Pk 

A .trl + fiP x ■fik 


K k 

P k-iVk 

if 

T k 

Pk- 1 X 

+ 'PkPk-i'Pk 

r l 

PU 

P k -xVkP k Pk-l 
\<rl + yjfcPfc-iVfc 

K k x fiP 


Pk 

1 ( P P k -i<PkV k P k -i \ 

* v +VlA-is?J 

1 (ft-. -*?-,) 

*9 


Table 1: The ('LSI and CLS2 implementation of the conventional RLS. 


In this table the (”) indicates that the same matrix- vector multiplication is performed 
for both implementations. 


2.3 The Square Root Covariance RLS Method (SCLS) 


This can be considered as a special case of the square- root -covariance Kalman-filter imple- 
mentation, also defined in Verhaegen and Van Dooren (1986). For the sytem matrices given 
in equation (15). the recursive scheme characterizing this implementation becomes 


/ 


V 


\/Al 


‘rk 


J S k . 1 


7T^-i 


Ux = 




0 


\ 


71 


s 


(16) 


* J 


pre-array 


post-array 


with the parameter update equation 

$k - dfr -1 -r — G k (y k - dl-\‘rlt) (IT) 

and Ui in equation (16) is an arbitrary orthogonal transformation that triangularizes the 


y/\ 

—rGk{y k - djLi^t-) 
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prearray. Here Sfe is the Choleski factor 3 of P^. defined as follows: 


Pk = Sk 5(. 


( 18 ) 


Remark 2: Another type of square root covariance RLS is the so-called UD-factorization 
algorithm (Biennan, 1976). The operational complexity of this scheme is slightly less than 
the presented SCLS scheme. However, from the point of view of numerical reliability, 
preliminary simulation tests performed in Verhaegen (1985) demonstrated a similar round- 
off error pattern between these two schemes. Therefore, the latter implementation will not 
be discussed further in this paper. □ 


2.4 The Square Root Information RLS Method (SILS) 


Similarly to the square-root-information Kalman-filter scheme, we can formulate for the LS 
problem recursions for the Choleski factor of the matrix Rk, defined in equation (4), and 
the parameter estimate d*, in a combined way. When the Choleski factor of Rk is defined 


as: 


this implementation becomes 

U2 



Rk = 

sATk-i 

v / A6- 

A 

<*k 

Mk 

(*k 


(19) 


\ 


( 20 ) 


Tk 6 

0 f, , 

with U 2 operating in a similar manner as U\. The update of the parameter estimates is 
given by: 

ih = K'Uk) (21) 


Remark 3: Other, so-called fast LS algorithms, such as the Levinson algorithm (Levinson* 
1947) are not included in this paper. This is because their numerical behavior has already 
been analyzed by an error perturbation study. For example for Levinson related algorithms 
this has been done in Cybenko. 1980. Furthermore, these algorithms impose a special 
structure upon the regression model, and this restriction is not considered in this paper. □ 
3 The '•Choleski factor" is often called the "square root.** In this paper the term "square root" is main- 
tained as far as the names of the square-root RLS schemes are concerned because of the familiarity that it 
has acquired. 
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3 Theoretical Error Analysis 


The numerical robustness with respect to round-off errors is analyzed for the four RLS 
schemes described in the previous section. Such an analysis can be split into three parts: 
(1) the round-off errors made in a single recursion, (2) the propagation of a single error in 
subsequent recursions and (3) the interaction and accumulation of the previous two error 
sources. In the paper of S. Ljung and L. Ljung (1985) only the second part has been treated. 
The following subsections present a much more refined analysis. 


3.1 Propagation of a Single Error Under oo Computational Precision 


In this subsection we consider the propagation of a single error at recursion instant k- 1 to 
subsequent recursions, assuming that no additional round-off errors are made. The following 
theorem describes this propagation for the CLSl and C'LS2 implementation. 

Theorem 1: Define the erroneous quantities, denoted by (.), at the k - 1 recursion instant 
as, 

P k—i = Pk-i + SPk-i and ih-i = d k-i + 6d k -\ (22) 

then these errors propagate to the next recursion instant k in the CLSl implementation as: 

6Pk = {(I ~ Kk'r[.)6Pk-i(I - Kky~',Y + \(6Pk-i - SP^^K', + 0(6 2 ) (23) 


64k = (/ - h'ky'k) 


/'V.f - : 


(A<rj; + ?' k P k - ifk) 


+ 0(1 r) 


and in the C'LS2 implementation as: 


1 


6Pk = j(I - Kk<Pk)6Pk-i(I - K\f' k y + 0(6 2 ) 


6h = (/ - Kky[.) 


bdk-\ t 


6Pk-i<rk 


( yk - r'k‘h-1 ) 


(\<ri + «ri.A--iy*r) 
where 0(6 2 ) indicates the order of magnitude of ij^P*._i|| 2 . 

Proof: (The proof will only be given for the CLSl implementation) 


+ 0 ( 6 2 ) 


(24) 

(25) 

(26) 


Consider the approximation. 


^ ff k + 'f'kPk-l'fk + 'r'k^Pk-l'fk A <t l + <fkPk 


2 ( 1 _ 'rk^Pk-X'fk \ , fit t2\ 

A 4 + tiPk-lfk ) 
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Then the error on K k updated in the CLSl becomes 

8K _ tPLvfk PLvPk /. _ ^k6Pk-v£k \ 0(6~) 

\<T k + if' k P k _i<p k \<T k + <p' k P k -xVk \ ^°k + PkPk-l'fk J 

Similar matrix-error analysis then results in the error propagation model (equation 23) for 

6P k . 

Following the same line of derivation, the error model (equation 24) for 6d k results. □ 

Corollary 1: The error propagation model (equation 23) demonstrates that for an ex- 
ponential forgetting factor A < 1, the effect of the loss of symmetry causes a blow up 
(divergence) in the error on P k . Therefore, when using the CLSl implementation it is abso- 
lutely necessary to maintain symmetry of the matrix P k to avoid divergence, no matter how 
accurately the computations are performed. This divergence phenomenon.' which is due to 
loss of symmetry, is well known, but not understood, in literature. Theorem 1, now gives a 
theoretical explanation of this phenomenon. This theorem furthermore indicates that such 
a loss of symmetry does not occur when using the CLS2 implementation. □ 

Corollary 2; Since for the moment we assume that no round-off errors are performed in 
the recursions, one inherently computes the same equations for the SCLS scheme as for the 
C'LS scheme. Therefore, starting with an error 6S k _i. which induces the following error on 

Pk- 1 : 

SPk. i = Sk-iSSi. x + + 6Sk.x6Si._x (27) 

and the error 6d k -\. the error propagation models ( equations 25-26) hold for the SCLS 
implementation. 

Remark 4: The same error models (equations 23-24) could be obtained by substituting 
the special system matrices (equation 15) in the error analysis results obtained for the 
conventional Kalman filter in Verhaegen and Van Dooren ( 1986). Repeating such an analysis 
particularly for the CLS. precisely indicates where the loss of symmetry comes from. □ 

For the SILS scheme, one inhtn nlly updates the information matrix R k according 
to equation (7). Here an error T k _ j, similarly represented as in equation (22), induces an 
error 6R k .\ given as: 

6Rk - 1 = T' k _x6T k - 1 + bTl.Jk-x + (28) 
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Again assuming no additional errors, this errors propagates according to: 

6R k = X 6Rk-i (29) 

Corollary 3: The error model (equation 29) immediately indicates that when using expo- 
nential forgetting and infinite precision , u information-matrix‘ , -type estimation schemes are 
exponentially stable. For A = 1 the error propagation does not increase. 

For the C’LSl, the CLS2 and the SCLS implementations, the matrix (I - Kk<p' k ) 
plays a crucial role in the propagation of both the errors 6 P k -\ and i • Therefore, let us 
inspect this matrix more closely. Substituting the expression for K k , obtained by comparing 
equations (6) and (8) and using equation (7) we find that: 

(l-K k <p' k ) = \R; l R k .-i (30) 

Now, consider the effect of an error induced at time instant ko < k on the computed 
quantities P k and When P k remains symmetric, this effect, becomes according to 
equation (25): 

6Pk = ^*(Mo)*JM(Mo)' + 0(* 2 ) (31) 

and according to equation (24): 

h-lto-l 

8d k = 4>{k.ko)8d ko + ^ <t>{ k, ko -I- j).pk 0 -ri+j + 0(b~) (32) 

j=o 

with /(, = — and <p(k,k Q +j) = U- =ko+j (I- A',v')- Assume now 

(A<r, + firi-iipi) 

that the conditions which guarantee that the residual (y, - ) becomes a zero-mean 

random variable hold. Then the second term in the right-hand side of equation ( 32) vanishes 
when we consider only the mean of as denoted by E( 6\h - )• Under this assumption, we 
will only analyze the propagation of E(8r) k ) in the sequel. 

The errors 6P ko and E(6d ko ) are attenuated if the transition matrix <p(k, k Q ) is a 
contraction. In the next theorem it will be shown that this is indeed the case if the regressor 
vector yk is persistently exciting. Let us first define this condition of the regressor vector. 

Definition 1: The regressor vector yv is persistently exciting over the observation interval 
k 0 < i < k when using an exponential forgetting factor A < 1, if the following condition is 
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fulfilled, 


al < £ W',* k ~' < PI 

r=fco 

for some positive constants a and 0. 


( 33 ) 


Theorem 2: When the regressor vector <p,- satisfies condition equation (33) over the interval 
(jfco, k), the transition matrix <f>(k, ko) is a contraction. 

Proof: 

The matrix R k is defined in equation (4) as, 

k 


Rk = Y. 

i=0 


Since k 0 < k we have that: 


R k = \ k - k *R ko + £ w\\ k ~' 

i~k o 

Therefore, when equation (33) holds, the matrix difference (R k - \ k '~ k ' 0 R ko ) is positive 
definite, or alternatively denoted: 

A k - k °R ko < R k 

When we now express <p(k.ko) as \ k ~ k ° R^ 1 R ko , using equation (30), <t>(k,k 0 ) is clearly 
positive definite. This completes the proof. 

Corollary 4: When the following conditions a re met. 

1. When the regressor vector is persistently exciting, and hence. R^ remains bounded 
(which can easily be shown). 

2. When A < 1. 

3. When no additional round-off errors are made in subsequent recursions. 


then t heorem 2 guarantees that, a single error 6P ko in equation (31) or E(St) kf> ) in equa- 
tion (32) decays exponentially when using the CLS2 and SCLS implementations. For the 
case A — 1. the contraction of R ^ 1 R ko guarantees that the error propagation remains stable. 
For the C'LSl implementation this holds when additionally the symmetry on P k is forced. 
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Another divergence phenomenon reported in literature is the loss of positive definite- 
ness of Pk . The question we might ask is how this phenomenon is related or interferes with 
the loss of symmetry discussed so far. Some insight into this complicated matter is gained 
by the following theorem. 

Theorem 3: When the errors on the parameter error covariance matrix Pk- 1 are symmetric 
and preserve its positive definiteness, the Pk updated by equation (10) remains positive 
definite. 


Proof: 

Since Pk-i is positive definite, it may be written as, 

Pk - 1 = XDX' 


(34) 


where XX' = /„ and D = diag(d \,- • • ,d„), with d,- > 0. Since X is of full rank, we might 
write the regressor vector pk as: 

Pk = Xp (35) 

with /i = (/i t , — , n„y. 

Substituting equations (34-35) in the update relationship for P we obtain: 


Pu = X{D - )A" 

Evaluating the matrix between brackets results in: 

(1 + n'Dp)Pk - 

<*1 ••• \ 


(36) 


JY 


-piPididi at 


\ -/<1 Hnd\d n -p 2 Undid" 


-pip„did„ 


ot, 


X' 


(37 


/ 


where o, = A <r|d,-+ V” =1 pjdjdj. The characteristic polynomial of the matrix between 
brackets in equation (37) can be evaluated for example by using t he symbolic manipulation 
software package by MACSYMA (1983) for different values of n. For all values of n, this 
polynomial can be written as: 


n-X 


(_l)n.n + ^(_i)i a ..» _ o with a, > 0 

i=0 

Obviously, this polynomial has only positive roots; therefore, Pk is positive definite. 


(38) 


II 



3.2 Round-Off Errors Made in a Single Recursion 


In the second part of the theoretical error analysis, we study the round-off errors made in 
a single recursion. The following theorem specifies bounds for the errors on the quantities 
updated each recursion in the four RLS schemes under investigation in this paper. 

Theorem 4: Denoting the norms of absolute errors caused by round off during the con- 
struction of P*, i?k, 5fc, Tfe, and P* by A p, A,j, As, At, and A p, respectively, we obtain 
the following upperbounds (where norms are 2-norms), 

1. CLS1 and CLS2 

A p < <i||Pfc|| 

A* < € 2 (||^|| + ||A' fc ||.||y fr ||) 

2. SCLS 

As < «3||S*||/c<W^l 

A P < €i\\Pk\\/ COS<i>\ 

A* < es(l|ty + IIAUM) 

3. S1LS 

At < ««||TV||/coa<p2 

A R < C;\\Rk\\lc03<i>2 

a p < <MPk)\micos<t> 2 

Aj < < 9 [K 2 (T k )\s k .\ + if(r*)Pfc|| + |f*|/co*fc] 
where «(.) denotes the condition number of the matrix (.), «■, are constants close to the 

machine precision t, cos<j>i are defined as follows: 

cos*, = II 11/ 11(67. | $*] I! 

cos<t> 2 = \\T k \\/ T '~ l 
I \ A- / ^ 

cos<t>3 = if*!/ 

which are usually close to one and s k is defined in equation (20). 

Proof: The proof would he a repetition of the one given in Verhaegen and Van Dooren 
(1986) for the corresponding Kalman filter implementations, taking into account the special 
system matrices (equation 15). Therefore, the reader is referred to this paper. 
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3.3 Accumulation and Interaction of the Round-Off Errors 


The third part of the analysis combines the local error upper bounds given in theorem 4 
with the propagation of a single error, to yield bounds on the total error for the different 
RLS schemes. The total error is denoted by the prefix Stot instead of 6. Here we make 
the assumption that the total error at some instant k is the sum of the propagation of the 
previous errors (when introducing no additional round-off errors), plus the local errors made 
during these recursions. 

For the CLS2, SCLS implementation the total error in the matrix i\. then satisfies: 

iiwwi < + s; (39) 

In this equation k m is the nearest time instant, smaller than k, for which 4>(k, k m ) is a 
contraction and is given as: 

5? = £ «<ll«ll (40) 

This represents the accumulation of the local' round-off errors made each recursion in the 
calculation of P, in the interval (k m ,k). When P; remains bounded, as imposed in condition 
1 of corollary 4, also remains bounded. 

In a similar way, the mean of the total error on the parameter estimates becomes: 

\\6tat*k\\ < |ld(*.* M )l!.||* ta ,d*JI + (41) 

Both equations (39) and (41) model the propagation of the total round-off errors in the 
CLS2 and SCLS implementation by a linear model with bounded input signal. 

Corollary 5: Error models (equations 39-41) indicate that when the regressor vector is 
persistently exciting, as defined in definition 1. the round-off errors made in the CLS2 
and SCLS implementations remain bounded. Furthermore, parts 1 and 2 of theorem 4 
show that the error models for both these implementations are identical. Therefore, only 
considering the numerical robustness no preference should be given to either of them. This 
latter condition generally does not hold for the CLSl implementation since the “ciires” to 
symmetrize P*. induce larger errors. This assertion has been confirmed in the experimental 
analysis made in Verhaegen and Van Dooren (1986) for Kalman filter implementations. 
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We now address the loss of positive definiteness under finite arithmetic precision. 
Using the assumption stated in this subsection to describe the total error on Pk , the following 
relationship results: 

(Pk + 6Pk) - AP < (P k + StotPk) < (Pk + SPk) + A P (42) 

with A P according to theorem 4 bounded as ||AP[j < e||(P* + SPk) ||. 

According to theorem 3, the matrix (Pk + SPk) remains positive definite and we can 
arrange its eigenvalues as: 

\x(Pk + SPk) > A t (P k + SP k .) > •■• > A n (P k + SPk) > 0 (43) 

The possibility of negative eigenvalues of (Pk'+ S tot Pk) arises with the lower bound 
in equation (42). Insight in the effect of the perturbation A P on the eigenvalues of (Pk + 
6P k ) - A P is given is given by the following lemma, taken from Wilkinson, pp. 101-102 
(1965). 

Lemma 1: If A and (A - E) are n by n symmetric matrices, then 

A n (.4 - E) > A„(A) - X\(E) 


Applying this lemma to the lower bound of equation (42). yields the following bound 
on the smallest eigenvalue of (Pk + S tot Pk): 

A „(Pk + StotPk) > A „(Pk + SP k .) - <A ,(/>* + SPk) (44) 

Therefore, this eigenvalue will be positive if 

1 ,45, 

This condition indicates that (Pk -r -SPk) is numerically not singular, as defined in Wilkinson 
(1965). 

Corpllary 6: If the recursion of the parameter covariance matrix Pk-i in equation (10) 
under finite precision preserves the symmetry of this matrix and this matrix is numerically 
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not singular, as stated by equation (45), than the updating via equation ( 10) can not destroy 
the positive definiteness of the parameter covariance matrix. 

From the computational scheme equation (20) of the SILS we clearly observe that 
it is not necessary to calculate the matrix P and the parameter estimate x) to continue 
recursions. Therefore, only the propagation of the total error on the information matrix Rk 
is of interest here. This model result from equation (29) and part 3 of theorem 4: 

||^tot.fifc|| < All^i^ll + e 7 ||Jfc-i|| (46) 

When the observations are persistently exciting, Rk increases; therefore, the only thing of 
interest here are the relative errors. If these are denoted by 6 T tot , equation (46) becomes 

\\r tot Rk\\<im 0 tRk-i\\ + <7 (47) 


where 7 = ||J*k-ill/Pfcll- 

Corollary 7s When- the observations are persistently exciting, the scalar 7 remains smaller 
than 1. Therefore, the linear relative error model (equation 47) shows that the computations 
with “information” type of RLS implementations will remain exponentially stable for A < 1. 
Each recursions, the quantities of interest Pk-i and may be computed from 7\._i and 
ik-i’ According to the bounds given in part 3 of theorem 4, these computations may be 
deteriorated by the condition number of Pk - 1 and Tk- 1. 

4 Experimental Evaluation 

The purpose of this section is to experimentally validate the corollaries made in the theo- 
retical error analysis of section 3. 

To restrict such a verification study, we refer to a similar study made in Verhae- 
gen and Van Dooren for Kalman filter implementations as an experimental verification of 
corollaries 2,3.5. and 7. 

Of particular interest to the scientist dealing with RLS problems is the numerical 
robustness of the conventional RLS scheme, since this scheme is very attractive in a broad 
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class of applications. In the present analysis corollaries 1,4, and 5 focus on this topic. More 
precisely, these three corollaries provide insight into the following questions: 

1. Whether a "wrong” way of implementing the conventional RLS has been the primary 
cause for its “high” sensitivity to the loss of symmetry (corollary 1)? 

2. What is the influence of persistent excitation on the numerical robustness of the RLS 
schemes (corollary 4)? 

3. What is the interference between the loss of symmetry and the loss of positive defi- 
niteness of Pit (corollary 6)? 

To experimentally evaluate these three corollaries, a mixed-precision simulation study 
is performed. Here the single-precision quantities, as denoted by (.). represent the erroneous 
quantities, and the double-precision quantities are assumed to be error free. In this exper- 
imental analysis we will only focus on the errors on P;. and with the above convention the 
total error on this quantity, as denoted by 6 tot in equation (39), is approximated by: 

11^^11 = 11^-^11 (48) 

The regression model taken as a test vehicle in this analysis has the following form: 

4 

Vk = (49) 

«=i 

Two different time histories of [yc- «<,?*•( 1),* • • .^(4)]. for k - 1. • • • .300. are used in the fol- 
lowing tests. Figure 1 displays these time histories, respectively, as measurement sequence 1 
and 2. 


The first, test evaluates corollary 1. For this purpose use is made of measurement 
sequence 1. which guarantees the presistentlv excitation condition to be fulfilled during 
the whole test. Furthermore, an exponential forgetting factor A = 0.95 was taken. The 
approximate error, equation (48), the loss of symmetry as computed by ||P<. - P^|| and 
||P;.j| are plotted for the CLSl implementation in Fig. 2a and for the C’LS2 implementation 
in Fig. 2b. These results clearly confirm corollary 1. 

In the second test corollary 6 is evaluated. The same experimental condition as in 
the previous test, are taken. Figures 3a and 3b, respectively, show the evolution of the 
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smallest eigenvalue and the third diagonal element of P k . Statistically these quantities 
have to be positive. In the error analysis, it has been reconfirmed that this remains true as 
long as P k remains symmetric and numerically not singular. The first condition is easily 
violated by the CLSl implementation; therefore, such a guarantee can not be stated for 
this implementation. This is confirmed in Fig. 3. Furthermore, such a loss of positive 
definiteness didn’t occur with the CLS2 implementation for all the test performed in this 
simulation study. From Fig. 2b, we clearly see that the loss of symmetry is of different 
orders of magnitude smaller than the errors on Pk- Therefore, the asymmetric part of AP 
in equation (42) can be neglected so that corollary 6 holds as well for the G'LS2 as for the 
SC'LS implementation. 

The evaluation of corollary 4 in the third test makes use of the measurement sequence 
2 in Fig. 1. From this figure we observe that the regressor vector pk is not persistently 
exciting during the time intervals (1,40) and (81,300). The impact of this on the round-off 
eiror propagation in the CLS2 implementation is pictured in Fig. 4. This figure clarifies 
that during the time intervals of lack of persistent excitation the error on P k and the loss of 
symmetry increase “linearly.” At the same time, ||Pj,|| increases also linearly. Furthermore, 
the rate of increase is exactly the same in both cases, so that relatively no loss of precision 
occurs. Therefore; the practical impact of the lack of persistent excitation causing the so- 
called “burst” phenomenon, see Hughes and Jacobs (1974). is far more important than the 
numerical behavior of the considered ELS schemes which preserve the symmetry of Pk- 


5 Concluding Remarks 

In this paper, the numerical robustness of four different RLS schemes is analyzed by means 
of a theoretical error analysis. Apart from reconfirming the insights obtained in a similar 
analysis for Kalman filter implementations, given in Verhaegen and Van Dooren (1986). 
These are. respectively: 

1. The conventional ELS (C'LS) which preserves symmetry of the parameter error co- 
variance matrix Pk yields the same accuracy as the SC'LS for all the corresponding 
quantities updated in both schemes. 
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2. The SILS is also numerically stable. The accuracy of computing Pk and the LS 
estimates is however penalized when the condition number of Pk is high. However, it 
is also not necessary to compute these quantities in order to continue the recursions 
with the SILS. Therefore, such large errors on Pk and the least squares estimates do 
not accumulate. 

The new insights, of particular interest for applications with RLS, of this analysis are: 

1. The CLS has been “ incorrectly " implemented, making it very to the so-called loss 
of symmetry phenomenon. Correct implementation of this scheme overcomes this 
deficiency. 

2. The persistency of excitation of the regressor vector in the regression model is required 
to guarantee the boundedness of the round-off errors. When this condition is violated, 
the errors on Pk will increase, however, with same rate of increase as ||P fr ||. Therefore, 
the increase of the errors is of minor importance. 

3. The loss of positive definiteness is caused by the loss of symmetry 
These conclusions have all been confirmed in an experimental verification study. 



References 


Astrom, K.J., U. Borisson, L. Ljung, and B. Wittenmark (1977). Theory and applications 
of self-tuning regulators. Automatica, 13 , pp. 457-476. 

Bierinan, G.J. (1976). Factorization methods for discrete sequential estimation. Academic 
Press , New York. 

Cybenko, G. ( 1980). The Numerical Stability of the Levison- Durbin algorithm for Toeplitz 
systems of equations. SIAM J. Sci. Stat. Comput.. 1, 3. pp. 303-319. 

Hagglund, T. (1983). New estimation techniques for adaptive control. Ph.D. dissertation , 
Lund Institute of Technology, Sweden. 

Hughes, D.J. and O.L. Jacobs (1974). Turn-off, escape and probing in non-linear stochastic 
control. IFAC Symp. on Stochastic control, pp. 343-352. 

Kulhavy, R. (1985). Restricted Exponential forgetting in Real-time identification. Seventh 
IFAC Symposium on Syst. Identification and Parameter Estimation , York. England. 

Levinson, N. (1947). The Weiner RMS (root mean square) error criterion in filter design 
and prediction. J. Math. Phys., 25 , pp. 261-278. 

Ljung, L. and T. Soderstrom (1983). Theory and practice of recursive identification. MIT 
Press. Massachusetts. 

Ljung, S. and L. Ljung (1985). Error propagation properties of recursive least squares 
adaptation algorithms. Automatica. 21. 2. pp. 157-167. 

MAC’SYMA Reference Manual (1983). The Mathlab Group Laboratory, Massachusetts 
Inst, of Tech., tenth ed. 

Verhaegen, M.H. and P. Van Dooren (1986). Numerical aspects of different Kalman filter 
implementations. IEEE Trans. Aut. C’ontr ., AC-31, Vol. 10, pp. 907-917, Oct. 1986. 

Verhaegen, M.H. (1985). A new class of algorithms in linear system theory, with appli- 
cation to real-time aerodynamic model identification. Ph.D. dissertation , Catholic 
University Leuven. Leuven, Belgium. 

Wilkinson. J.H. (1965). The Algebraic Eigenvalue Problem. Clarendon Press. Oxford. 


19 



Figure Captions 


Fig. 1 : Two different measurement sequences for the variables (j/fc, y? u ) in the regression 
model (1). 

Fig. 2~. Propagation of round-off errors in two implementations of the conventional RLS 
scheme (A = 0.95, measurement sequence 1). 

Fig. 3 : Evolution of certain eigenvalues and diagonal elements of Pk in the C'LS 1 imple- 
mentation (A = 0.95, measurement sequence 1). 

Fig. 4: Influence off the persistency of excitation of the regressor vector y? ( h) on the round- 
off propagation in CLS 2 implementation (A = 0.95, measurement sequence 2). 
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(b) THE CLS 2 IMPLEMENTATION 

Fig. 2 
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Fig. 4 


24 


NASA 

Nation, it Aeronautics and 
Space Admintsl ration 


Report Documentation Page 


1. Report No. 

NASA TM- 100055 


2. Government Accession No. 


3. Recipient's Catalog No. 


4. Title and Subtitle 

Round-off Error Propagation in Four Generally Applicable, 
Recursive, Least-Squares-Estimation Schemes 


5. Report Date 

December 1987 


6. Performing Organization Code 


7. Author(s) 

M. H. Verhaegen 


8. Performing Organization Report No. 

A-88050 


9. Performing Organization Name and Address 

Ames Research Center 
Moffett Field, CA 94035 


10. Work Unit No. 

505-66-11 


11. Contract or Grant No. 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


13. Type of Report and Period Covered 

Technical Memorandum 


14. Sponsoring Agency Code 


15. Supplementary Notes 

Point of Contact: Ralph Bach, Jr., Ames Research Center, MS 210-9, Moffett Field, CA 94035 
(415) 694-5429 or FTS 464-5429 


16. Abstract 


In this paper the numerical robustness of four generally applicable, recursive, least-squares- 
estimation schemes is analyzed by means of a theoretical round-off propagation study. This study 
highlights a number of practical, interesting insights of widely used recursive least-squares schemes. 
These insights have been confirmed in an experimental verification study as well. 


17. Key Words (Suggested by Author(s)) 

Recursive least squares 
Divergence phenomenon 
Theoretical error analysis 


18. Distribution Statement 

Unclassified - Unlimited 

Subject Category - 64 


19. Security Classif. (of this report) 

20. Security Classif. (of this page) 

21. No. of pages 

22. Price 

Unclassified 

Unclassified 

25 

A02 


NASA FORM 1626 OCT 86 


For sale by the National Technical Information Service, Springfield, Virginia 22161 





